﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using System.Collections.Specialized;
using System.Text.RegularExpressions;

namespace CS281B_Project
{
    class Program
    {
        static void Main(string[] args)
        {
            #region Argument
            int K = 10;
            int seedNum = 4;
            string confName = "config.txt";
            string modelName = "baseline";
            bool willOutputMean = false;
            bool willPredict = true;
            bool willTrain = true;
            string importModel = null;
            Arguments options = new Arguments(args);
            if (options["K"] != null)
                K = int.Parse(options["K"]);
            if (options["ConfName"] != null)
                confName = options["ConfName"];
            if (options["ModelName"] != null)
                modelName = options["ModelName"];
            if (options["SeedNum"] != null)
                seedNum = int.Parse(options["SeedNum"]);
            if (options["DoNotTrain"] != null)
                willTrain = false;
            if (options["DoNotPredict"] != null)
                willPredict = false;
            if (options["OutputMean"] != null)
                willOutputMean = true;
            if (options["ImportModel"] != null)
                importModel = options["ImportModel"];

            #endregion

            Console.WriteLine(string.Format("Training begins. K={0}, conf={1}, model={2}, seedNum={3}, outputMean={4}, willTrain={5}, willPredict={6}, importModel={7}", K, confName, modelName, seedNum, willOutputMean, willTrain, willPredict, importModel == null ? "None" : importModel));
            Trainer trainer = new Trainer(K, seedNum);
            trainer.ReadConfigure(confName);
            trainer.Initialize();

            int[] index = new int[150];
            int[] trainSet = new int[120];
            int[] testSet = new int[30];

            Random rnd = new Random(1);

            for (int i = 1; i <= 150; i++)
            {
                if (i % 5 != 0)
                {
                    trainSet[i - i / 5 - 1] = i;
                }
                else
                {
                    testSet[i / 5 - 1] = i;
                }
            }
            shuffle(trainSet, rnd);
            shuffle(testSet, rnd);

            if (importModel != null)
                trainer.ImportModel(importModel); // (@"models/model_K0_train.txt");

            if (willTrain)
            {
                for (int iter = 0; iter < 1; iter++)
                {
                    for (int i = 0; i < 15; i++)
                    {
                        Console.WriteLine("Group {0}:", i + 1);
                        for (int j = 0; j < 8; j++)
                        {
                            trainer.Train(trainSet[i * 8 + j]);
                        }

                        trainer.ExportModel(@"models/model_" + modelName + ".txt");
                        //trainer.Test(testSet[i * 2]);
                        //trainer.Test(testSet[i * 2 + 1]);

                    }
                }

                for (int i = 0; i < 30; i++)
                {
                    trainer.Train(testSet[i]);
                }

                trainer.ExportModel(@"models/model_" + modelName + ".txt");
            }

            if (willPredict)
                trainer.Predict(151, @"models/prediction_" + modelName + ".csv", 0);

            if (willOutputMean)
                for (int i = 1; i <= 151; i++)
                {
                    trainer.OutputMean(i, modelName);
                }

            //Console.ReadKey();
        }

        public static void shuffle(int[] set, Random seed)
        {
            for (int i = set.Length - 1; i >= 1; i--)
            {
                int r = seed.Next(i);
                int t = set[i];
                set[i] = set[r];
                set[r] = t;
            }
        }
    }

    public class Trainer
    {
        int K = 0;
        double var_uv = 1.0;
        double var_w = 2.0;
        int AMPIterNum = 3;
        int featureNum;
        //double time_measure = 1;
        Feature dim = new Feature();
        Feature index = new Feature();
        Feature group = new Feature();

        Gaussian[,][] u;
        Gaussian[,][] v;
        Gaussian[][] w;
        Dictionary<int, Gaussian[][]> subweights = new Dictionary<int, Gaussian[][]>();

        double[] ll = new double[4];
        double[] impr = new double[4];

        Random seed;

        public Trainer(int K, int seedNum)
        {
            seed = new Random(seedNum);
            this.K = K;
        }
        public void ImportPrior(String path)
        {
            StreamReader sr = new StreamReader(path);
            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.x.Length; j++)
                {
                    for (int k = 0; k < dim.x[j]; k++)
                    {
                        String[] items = sr.ReadLine().Split(new string[] { "(", ")", "," }, StringSplitOptions.RemoveEmptyEntries);
                        if (items.Length == 2)
                        {
                            u[i, j][k] = new Gaussian(double.Parse(items[0]), var_uv / group.x.Length / K);
                        }
                    }
                }
            }

            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.y.Length; j++)
                {
                    for (int k = 0; k < dim.y[j]; k++)
                    {
                        String[] items = sr.ReadLine().Split(new string[] { "(", ")", "," }, StringSplitOptions.RemoveEmptyEntries);
                        if (items.Length == 2)
                        {
                            v[i, j][k] = new Gaussian(double.Parse(items[0]), var_uv / group.y.Length / K);
                        }
                    }
                }
            }

            for (int j = 0; j < dim.p.Length; j++)
            {
                for (int k = 0; k < dim.p[j]; k++)
                {
                    String[] items = sr.ReadLine().Split(new string[] { "(", ")", "," }, StringSplitOptions.RemoveEmptyEntries);
                    if (items.Length == 2)
                    {
                        w[j][k] = new Gaussian(double.Parse(items[0]), var_w / group.p.Length);
                    }
                }
            }
            sr.Close();
        }

        public void ImportModel(String path)
        {
            StreamReader sr = new StreamReader(path);
            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.x.Length; j++)
                {
                    for (int k = 0; k < dim.x[j]; k++)
                    {
                        String[] items = sr.ReadLine().Split(new string[] { "(", ")", "," }, StringSplitOptions.RemoveEmptyEntries);
                        if (items.Length == 2)
                        {
                            u[i, j][k] = new Gaussian(double.Parse(items[0]), double.Parse(items[1]));
                        }
                    }
                }
            }

            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.y.Length; j++)
                {
                    for (int k = 0; k < dim.y[j]; k++)
                    {
                        String[] items = sr.ReadLine().Split(new string[] { "(", ")", "," }, StringSplitOptions.RemoveEmptyEntries);
                        if (items.Length == 2)
                        {
                            v[i, j][k] = new Gaussian(double.Parse(items[0]), double.Parse(items[1]));
                        }
                    }
                }
            }

            for (int j = 0; j < dim.p.Length; j++)
            {
                for (int k = 0; k < dim.p[j]; k++)
                {
                    String[] items = sr.ReadLine().Split(new string[] { "(", ")", "," }, StringSplitOptions.RemoveEmptyEntries);
                    if (items.Length == 2)
                    {
                        w[j][k] = new Gaussian(double.Parse(items[0]), double.Parse(items[1]));
                    }
                }
            }
            sr.Close();
        }

        public void ExportModel(String path)
        {
            StreamWriter sw = new StreamWriter(path);
            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.x.Length; j++)
                {
                    for (int k = 0; k < dim.x[j]; k++)
                    {
                        if (u[i, j][k] != null)
                            sw.WriteLine("({0:G8},{1:G8})", u[i, j][k].mean, u[i, j][k].variance);
                        else
                            sw.WriteLine("()");
                    }
                }
            }

            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.y.Length; j++)
                {
                    for (int k = 0; k < dim.y[j]; k++)
                    {
                        if (v[i, j][k] != null)
                            sw.WriteLine("({0:G8},{1:G8})", v[i, j][k].mean, v[i, j][k].variance);
                        else
                            sw.WriteLine("()");
                    }
                }
            }

            for (int j = 0; j < dim.p.Length; j++)
            {
                for (int k = 0; k < dim.p[j]; k++)
                {
                    if (w[j][k] != null)
                        sw.WriteLine("({0:G8},{1:G8})", w[j][k].mean, w[j][k].variance);
                    else
                        sw.WriteLine("()");
                }
            }
            sw.Close();
        }

        public void Train(int N)
        {
            FastBinaryReader trainStream = new FastBinaryReader(@"features/feature_" + N + ".bin");
            FastBinaryReader biasStream = new FastBinaryReader(@"bias/mean_" + N + ".bin");
            Console.WriteLine("Training: File {0}", N);
            Feature current;
            while ((current = ReadNextLine(trainStream, biasStream)) != null)
            {
                for (int i = 0; i < current.impr; i++)
                {
                    double[] a = new double[2];
                    if (i < current.click)
                    {
                        a[0] = 0;
                        a[1] = 1;
                    }
                    else
                    {
                        a[0] = 1;
                        a[1] = -1;
                    }

                    UpdateWeights(current, a, seed);
                }
            }
            trainStream.Close();
            biasStream.Close();
        }

        public void Test(int N)
        {
            FastBinaryReader testStream = new FastBinaryReader(@"features/feature_" + N + ".bin");
            FastBinaryReader biasStream = new FastBinaryReader(@"bias/mean_" + N + ".bin");
            Console.WriteLine("Testing: File {0}", N);
            Feature current;
            while ((current = ReadNextLine(testStream, biasStream)) != null)
            {
                double p = returnCTR(current);
                impr[0] += current.impr;
                impr[current.position + 1] += current.impr;
                ll[0] += current.click * Math.Log(p) + (current.impr - current.click) * Math.Log(1 - p);
                ll[current.position + 1] += current.click * Math.Log(p) + (current.impr - current.click) * Math.Log(1 - p);
            }
            testStream.Close();
            biasStream.Close();
            Console.WriteLine("Log-Likelihood: {0}\t{1}\t{2}\t{3}", -ll[0] / impr[0], -ll[1] / impr[1], -ll[2] / impr[2], -ll[3] / impr[3]);
        }

        public void OutputMean(int N)
        {
            OutputMean(N, null);
        }

        public void OutputMean(int N, string modelName)
        {
            if (modelName == null)
                modelName = "new";
            if (!System.IO.Directory.Exists("bias_" + modelName))
                Directory.CreateDirectory("bias_" + modelName);
            FastBinaryReader input = new FastBinaryReader(@"features/feature_" + N + ".bin");
            FastBinaryReader biasStream = new FastBinaryReader(@"bias/mean_" + N + ".bin");
            StreamWriter output = new StreamWriter(@"bias_" + modelName + @"/mean_" + N + ".txt");
            Console.WriteLine("Outputing Mean: File {0}", N);

            Feature current;
            int counter = 0;
            double[] a = new double[2];
            a[0] = 1;
            a[1] = 0;

            while ((current = ReadNextLine(input, biasStream)) != null)
            {
                counter++;
                if (counter % 1000000 == 0)
                    Console.WriteLine(counter);

                Gaussian[] s = new Gaussian[K];
                Gaussian[] t = new Gaussian[K];
                Gaussian b = null;
                PreparePriors(current, s, t, ref b, null);
                Gaussian posterior = ApproximateMessagePassing(s, t, ref b, a);
                output.WriteLine("{0:N6}", posterior.mean);
            }
            input.Close();
            biasStream.Close();
            output.Close();
        }

        public void Predict(int N, String path, double t)
        {
            FastBinaryReader input = new FastBinaryReader(@"features/feature_" + N + ".bin");
            FastBinaryReader biasStream = new FastBinaryReader(@"bias/mean_" + N + ".bin");
            StreamWriter output = new StreamWriter(path);
            Console.WriteLine("Predicting: File {0}", N);

            Feature current;
            int counter = 0;
            int cut = 0;
            int[] impr = new int[3];
            double[] ll = new double[3];
            while ((current = ReadNextLine(input, biasStream)) != null)
            {
                counter++;
                if (counter % 1000000 == 0)
                    Console.WriteLine(counter);
                double p = returnCTR(current);
                if (p < t)
                {
                    p = 0;
                    cut++;
                }
                output.WriteLine("{0:N6}", p);
                output.Flush();
            }
            input.Close();
            biasStream.Close();
            output.Close();
        }

        public void ReadConfigure(String path)
        {
            StreamReader sr = new StreamReader(path);
            String line = sr.ReadLine();
            
            string[] items = line.Split(new string[] { " ", "\t" }, StringSplitOptions.RemoveEmptyEntries);
            string[] subitems;
            string[] subitems2;

            subitems = items[0].Split(new string[] { "[", "]", "," }, StringSplitOptions.RemoveEmptyEntries);
            subitems2 = items[0].Split(new string[] { "[", "]", ",", "|", ":" }, StringSplitOptions.RemoveEmptyEntries);
            dim.x = new int[subitems.Length];
            index.x = new int[subitems2.Length - subitems.Length];
            group.x = new int[subitems2.Length - subitems.Length];

            int dindex = 0;
            int gindex = 0;
            for (int i = 0; i < subitems.Length; i++)
            {
                string[] subsubitems = subitems[i].Split(new string[] { ":", "|" }, StringSplitOptions.RemoveEmptyEntries);
                for (int j = 0; j < subsubitems.Length - 1; j++)
                {
                    index.x[dindex] = int.Parse(subsubitems[j]);
                    group.x[dindex] = gindex;
                    dindex++;
                }
                dim.x[gindex] = int.Parse(subsubitems[subsubitems.Length - 1]);
                gindex++;
            }

            subitems = items[1].Split(new string[] { "[", "]", "," }, StringSplitOptions.RemoveEmptyEntries);
            subitems2 = items[1].Split(new string[] { "[", "]", ",", "|", ":" }, StringSplitOptions.RemoveEmptyEntries);
            dim.y = new int[subitems.Length];
            index.y = new int[subitems2.Length - subitems.Length];
            group.y = new int[subitems2.Length - subitems.Length];

            dindex = 0;
            gindex = 0;
            for (int i = 0; i < subitems.Length; i++)
            {
                string[] subsubitems = subitems[i].Split(new string[] { ":", "|" }, StringSplitOptions.RemoveEmptyEntries);
                for (int j = 0; j < subsubitems.Length - 1; j++)
                {
                    index.y[dindex] = int.Parse(subsubitems[j]);
                    group.y[dindex] = gindex;
                    dindex++;
                }
                dim.y[gindex] = int.Parse(subsubitems[subsubitems.Length - 1]);
                gindex++;
            }

            subitems = items[2].Split(new string[] { "[", "]", "," }, StringSplitOptions.RemoveEmptyEntries);
            subitems2 = items[2].Split(new string[] { "[", "]", ",", "|", ":" }, StringSplitOptions.RemoveEmptyEntries);
            dim.p = new int[subitems.Length];
            index.p = new int[subitems2.Length - subitems.Length];
            group.p = new int[subitems2.Length - subitems.Length];

            dindex = 0;
            gindex = 0;
            for (int i = 0; i < subitems.Length; i++)
            {
                string[] subsubitems = subitems[i].Split(new string[] { ":", "|" }, StringSplitOptions.RemoveEmptyEntries);
                for (int j = 0; j < subsubitems.Length - 1; j++)
                {
                    index.p[dindex] = int.Parse(subsubitems[j]);
                    group.p[dindex] = gindex;
                    dindex++;
                }
                dim.p[gindex] = int.Parse(subsubitems[subsubitems.Length - 1]) + 1;
                gindex++;
            }

            line = sr.ReadLine();
            featureNum = int.Parse(line);
            sr.Close();
        }

        public Feature ReadNextLine(FastBinaryReader sr, FastBinaryReader sr_bias)
        {
            Feature f = null;
            int[] items = new int[featureNum];
            if (!sr.EndOfStream())
            {
                sr.NextIntArray(featureNum, ref items);

                f = new Feature();
                f.click = items[0];
                f.impr = items[1];
                f.position = items[6];

                f.x = new int[group.x.Length];
                for (int i = 0; i < group.x.Length; i++)
                {
                    if (items[index.x[i]] == 0)
                        f.x[i] = -1;
                    else
                        f.x[i] = items[index.x[i]] % (dim.x[group.x[i]] - 1);
                }

                f.y = new int[group.y.Length];
                for (int i = 0; i < group.y.Length; i++)
                {
                    if (items[index.y[i]] == 0)
                        f.y[i] = -1;
                    else
                        f.y[i] = items[index.y[i]] % (dim.y[group.y[i]] - 1);
                }

                f.p = new int[group.p.Length];
                for (int i = 0; i < group.p.Length; i++)
                {
                    f.p[i] = items[index.p[i]] % (dim.p[group.p[i]] - 1);
                }

                f.q = sr_bias.NextDouble();
            }
            return f;
        }

        public Feature ReadNextLine(StreamReader sr, StreamReader sr_bias)
        {
            String line;
            Feature f = null;
            if ((line = sr.ReadLine()) != null)
            {
                string[] items = line.Split(new string[] { " ", "\t", "[", "]", "," }, StringSplitOptions.RemoveEmptyEntries);

                f = new Feature();
                f.click = int.Parse(items[0]);
                f.impr = int.Parse(items[1]);
                f.position = int.Parse(items[6]);

                f.x = new int[group.x.Length];
                for (int i = 0; i < group.x.Length; i++)
                {
                    if (int.Parse(items[index.x[i]]) == 0)
                        f.x[i] = -1;
                    else
                        f.x[i] = int.Parse(items[index.x[i]]) % (dim.x[group.x[i]] - 1);
                }

                f.y = new int[group.y.Length];
                for (int i = 0; i < group.y.Length; i++)
                {
                    if (int.Parse(items[index.y[i]]) == 0)
                        f.y[i] = -1;
                    else
                        f.y[i] = int.Parse(items[index.y[i]]) % (dim.y[group.y[i]] - 1);
                }

                f.p = new int[group.p.Length];
                for (int i = 0; i < group.p.Length; i++)
                {
                    f.p[i] = int.Parse(items[index.p[i]]) % (dim.p[group.p[i]] - 1);
                }

                f.q = double.Parse(sr_bias.ReadLine());
            }
            return f;
        }

        public void Initialize()
        {
            u = new Gaussian[K, dim.x.Length][];
            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.x.Length; j++)
                {
                    u[i, j] = new Gaussian[dim.x[j]];
                }
            }

            v = new Gaussian[K, dim.y.Length][];
            for (int i = 0; i < K; i++)
            {
                for (int j = 0; j < dim.y.Length; j++)
                {
                    v[i, j] = new Gaussian[dim.y[j]];
                }
            }

            w = new Gaussian[dim.p.Length][];
            for (int j = 0; j < dim.p.Length; j++)
            {
                w[j] = new Gaussian[dim.p[j]];
                w[j][dim.p[j] - 1] = new Gaussian(0, var_w / group.p.Length);
            }
        }


        public double returnCTR(Feature feature)
        {
            Gaussian[] s = new Gaussian[K];
            Gaussian[] t = new Gaussian[K];
            Gaussian b = null;
            double[] a = new double[2];
            a[0] = 1;
            a[1] = 0;

            PreparePriors(feature, s, t, ref b, null);
            Gaussian posterior = ApproximateMessagePassing(s, t, ref b, a);
            double p = GaussianIntegral.NormCdf(posterior.mean / Math.Sqrt(1));
            return p;
        }

        public void UpdateWeights(Feature feature, double[] a, Random seed)
        {
            Gaussian[] s = new Gaussian[K];
            Gaussian[] t = new Gaussian[K];
            Gaussian b = null;
            PreparePriors(feature, s, t, ref b, seed);//Generate Priors: s,t,b

            Gaussian[] s_prior = new Gaussian[K];
            Gaussian[] t_prior = new Gaussian[K];
            Gaussian b_prior;
            for (int i = 0; i < K; i++)
            {
                s_prior[i] = new Gaussian(s[i].mean, s[i].variance);
                t_prior[i] = new Gaussian(t[i].mean, t[i].variance);
            }
            b_prior = new Gaussian(b.mean, b.variance);

            ApproximateMessagePassing(s, t, ref b, a);//Message Passing
            Gaussian b_likelihood = Gaussian.divide(b, b_prior);


            for (int j = 0; j < feature.p.Length; j++)//Update Weights: w
            {
                if (feature.p[j] == -1)
                    continue;

                int g = group.p[j];
                Gaussian w_likelihood = new Gaussian(b_likelihood.mean - b_prior.mean + w[g][feature.p[j]].mean,
                        b_likelihood.variance + b_prior.variance - w[g][feature.p[j]].variance);

                w[g][feature.p[j]] = Gaussian.multiply(w[g][feature.p[j]], w_likelihood);
                w[g][dim.p[g] - 1] = Gaussian.multiply(w[g][dim.p[g] - 1], w_likelihood);
                //w[j][feature.p[j]].variance *= time_measure;
                //w[j][dim.p[j] - 1].variance *= time_measure;
            }

            for (int i = 0; i < K; i++)//Update Weights: u,v
            {
                Gaussian s_likelihood = Gaussian.divide(s[i], s_prior[i]);
                for (int j = 0; j < feature.x.Length; j++)
                {
                    if (feature.x[j] == -1)
                        continue;

                    int g = group.x[j];
                    Gaussian u_likelihood = new Gaussian(s_likelihood.mean - s_prior[i].mean + u[i, g][feature.x[j]].mean,
                        s_likelihood.variance + s_prior[i].variance - u[i, g][feature.x[j]].variance);
                    u[i, g][feature.x[j]] = Gaussian.multiply(u[i, g][feature.x[j]], u_likelihood);
                    //u[i, g][feature.x[j]].variance *= time_measure;
                }

                Gaussian t_likelihood = Gaussian.divide(t[i], t_prior[i]);
                for (int j = 0; j < feature.y.Length; j++)
                {
                    if (feature.y[j] == -1)
                        continue;

                    int g = group.y[j];
                    Gaussian v_likelihood = new Gaussian(t_likelihood.mean - t_prior[i].mean + v[i, g][feature.y[j]].mean,
                        t_likelihood.variance + t_prior[i].variance - v[i, g][feature.y[j]].variance);
                    v[i, g][feature.y[j]] = Gaussian.multiply(v[i, g][feature.y[j]], v_likelihood);
                    //v[i, g][feature.y[j]].variance *= time_measure;
                }
            }
        }

        public void PreparePriors(Feature feature, Gaussian[] s, Gaussian[] t, ref Gaussian b, Random seed)
        {

            double meansum;
            double varsum;
            for (int i = 0; i < K; i++)//Calculate Prior
            {
                meansum = 0;
                varsum = 0;
                for (int j = 0; j < feature.x.Length; j++)
                {
                    if (feature.x[j] == -1)
                        continue;
                    int g = group.x[j];
                    if (u[i, g][feature.x[j]] == null)
                        u[i, g][feature.x[j]] = new Gaussian((NextDouble(seed) - 0.5) * 2 / Math.Sqrt(feature.x.Length * K), var_uv / feature.x.Length / K);
                    meansum += u[i, g][feature.x[j]].mean;
                    varsum += u[i, g][feature.x[j]].variance;
                }
                s[i] = new Gaussian(meansum, varsum);

                meansum = 0;
                varsum = 0;
                for (int j = 0; j < feature.y.Length; j++)
                {
                    if (feature.y[j] == -1)
                        continue;
                    int g = group.y[j];
                    if (v[i, g][feature.y[j]] == null)
                        v[i, g][feature.y[j]] = new Gaussian((NextDouble(seed) - 0.5) * 2 / Math.Sqrt(feature.y.Length * K), var_uv / feature.y.Length / K);
                    meansum += v[i, g][feature.y[j]].mean;
                    varsum += v[i, g][feature.y[j]].variance;
                }
                t[i] = new Gaussian(meansum, varsum);
            }
            meansum = 0;
            varsum = 0;
            for (int j = 0; j < feature.p.Length; j++)
            {
                if (feature.p[j] == -1)
                    continue;
                int g = group.p[j];
                if (w[g][feature.p[j]] == null)
                    w[g][feature.p[j]] = new Gaussian(w[g][dim.p[g] - 1].mean, var_w / feature.p.Length);
                else if (w[g][feature.p[j]].variance >= var_w / feature.p.Length - 1e-8)
                    w[g][feature.p[j]].mean = w[g][dim.p[g] - 1].mean;
                meansum += w[g][feature.p[j]].mean;
                varsum += w[g][feature.p[j]].variance;
            }
            b = new Gaussian(meansum + feature.q, varsum);
        }

        public Gaussian ApproximateMessagePassing(Gaussian[] s, Gaussian[] t, ref Gaussian b, double[] a)
        {
            Gaussian[] s_prior = new Gaussian[K];
            Gaussian[] t_prior = new Gaussian[K];
            Gaussian[] star_to_s = new Gaussian[K];
            Gaussian[] star_to_t = new Gaussian[K];
            Gaussian[] star_to_z = new Gaussian[K];
            Gaussian[] z_to_star = new Gaussian[K];

            for (int i = 0; i < K; i++)
            {
                s_prior[i] = new Gaussian(s[i].mean, s[i].variance);
                t_prior[i] = new Gaussian(t[i].mean, t[i].variance);
            }

            double meansum = 0;
            double varsum = 0;
            Gaussian likelihood = null;
            Gaussian r_posterior = null;
            for (int iter = 0; iter < AMPIterNum; iter++)
            {
                meansum = 0;
                varsum = 0;
                double[] st = new double[K];
                double[] s2 = new double[K];
                double[] t2 = new double[K];

                for (int i = 0; i < K; i++)
                {
                    st[i] = s[i].mean * t[i].mean;
                    s2[i] = s[i].variance + s[i].mean * s[i].mean;
                    t2[i] = t[i].variance + t[i].mean * t[i].mean;
                    star_to_z[i] = new Gaussian(st[i], s2[i] * t2[i] - st[i] * st[i]);
                    meansum += star_to_z[i].mean;
                    varsum += star_to_z[i].variance;
                }
                Gaussian r_prior = new Gaussian(meansum + b.mean, varsum + b.variance);
                r_posterior = new Gaussian(r_prior.mean, r_prior.variance);
                r_posterior.update(a);
                likelihood = Gaussian.divide(r_posterior, r_prior);
                for (int i = 0; i < K; i++)
                {
                    z_to_star[i] = new Gaussian(likelihood.mean - b.mean - meansum + star_to_z[i].mean, likelihood.variance + b.variance + varsum - star_to_z[i].variance);
                    star_to_s[i] = new Gaussian(z_to_star[i].mean * t[i].mean / t2[i], z_to_star[i].variance / t2[i]);
                    star_to_t[i] = new Gaussian(z_to_star[i].mean * s[i].mean / s2[i], z_to_star[i].variance / s2[i]);
                    s[i] = Gaussian.multiply(s_prior[i], star_to_s[i]);
                    t[i] = Gaussian.multiply(t_prior[i], star_to_t[i]);
                }
            }
            Gaussian b_likelihood = new Gaussian(likelihood.mean - meansum, likelihood.variance + varsum);
            b = Gaussian.multiply(b, b_likelihood);

            return r_posterior;
        }

        public static double NextDouble(Random seed)
        {
            if (seed != null)
                return seed.NextDouble();
            else
                return 0;
        }
    }

    public class FastBinaryReader
    {
        const int BufferSize = 1 << 20;
        byte[] buffer;
        int bufferContentLength;
        int curIndex;
        FileStream fileStream;

        public FastBinaryReader(string filePath)
        {
            fileStream = new FileStream(filePath, FileMode.Open, FileAccess.Read, FileShare.Read, BufferSize);
            buffer = new byte[BufferSize];
            curIndex = 0;
            bufferContentLength = 0;
        }

        public double NextDouble()
        {
            if (curIndex == bufferContentLength)
            {
                bufferContentLength = fileStream.Read(buffer, 0, BufferSize);
                curIndex = 0;
            }
            curIndex += 8;
            return BitConverter.ToDouble(buffer, curIndex - 8);
        }

        public int NextInt()
        {
            if (curIndex == bufferContentLength)
            {
                bufferContentLength = fileStream.Read(buffer, 0, BufferSize);
                curIndex = 0;
            }
            curIndex += 4;
            return BitConverter.ToInt32(buffer, curIndex - 4);
        }

        public void NextIntArray(int count, ref int[] result)
        {
            for (int i = 0; i < count; i++) 
            {
                result[i] = NextInt();
            }
        }

        public bool EndOfStream()
        {
            if (curIndex == bufferContentLength)
            {
                bufferContentLength = fileStream.Read(buffer, 0, BufferSize);
                curIndex = 0;
                return (bufferContentLength == 0);
            }
            else
            {
                return false;
            }
        }

        public void Close()
        {
            fileStream.Close();
        }
    }


    public class Feature
    {
        //public int adID;
        public int click;
        public int impr;
        public int position;
        public int[] x;
        public int[] y;
        public int[] p;
        public double q = 0;
    }

    public class Gaussian
    {
        public double mean = 0;
        public double variance = 1;

        public Gaussian(double x1, double x2)
        {
            mean = x1;
            variance = x2;
        }

        public static Gaussian multiply(Gaussian a, Gaussian b)
        {
            return new Gaussian((a.mean * b.variance + b.mean * a.variance) / (a.variance + b.variance),
                a.variance * b.variance / (a.variance + b.variance));
        }

        public static Gaussian divide(Gaussian a, Gaussian b)
        {
            if (a.variance >= b.variance)
                return new Gaussian(0, 1);
            else
                return new Gaussian((a.mean * b.variance - b.mean * a.variance) / (b.variance - a.variance),
               a.variance * b.variance / (b.variance - a.variance));
        }

        public static double[] Integral(double mu, double delta, int k)
        {
            double[] moments = new double[3];
            double muZX, deltaZX, deltaXZ, dk, ek, vk, s;

            if (k == 0)
            {
                moments[0] = 1;
                moments[1] = mu;
                moments[2] = delta + mu * mu;
                return moments;
            }

            deltaXZ = delta + 1;       //(7)
            dk = GaussianIntegral.NormCdf(mu / Math.Sqrt(deltaXZ));       //(8)
            ek = mu + Math.Sqrt(deltaXZ) / Math.Sqrt(2 * Math.PI) / dk * Math.Exp(-mu * mu / 2 / deltaXZ);      //(9)
            vk = deltaXZ + mu * ek - ek * ek;   //(10)
            muZX = (deltaXZ * ek - vk * mu) / (deltaXZ - vk);      //(11)
            deltaZX = vk * deltaXZ / (deltaXZ - vk) + 1;   //(12)

            s = dk * Math.Sqrt(2 * Math.PI * (deltaXZ + deltaZX - 1)) * Math.Exp(Math.Pow(mu - muZX, 2) / 2 / (deltaXZ + deltaZX - 1));
            moments[0] = s * Math.Exp(-Math.Pow(mu - muZX, 2) / 2 / (deltaZX + delta)) / Math.Sqrt((2 * Math.PI) * (deltaZX + delta));
            moments[1] = moments[0] * (delta * muZX + deltaZX * mu) / (delta + deltaZX);
            moments[2] = moments[0] * delta * deltaZX / (delta + deltaZX) + moments[1] * moments[1] / moments[0];

            return moments;
        }

        public void update(double[] a)
        {
            int length = a.Length;
            double auSum = 0, acSum = 0, avSum = 0;
            double[] moments;
            for (int k = 0; k < a.Length; k++)
            {
                moments = Integral(mean, variance, k);
                acSum += a[k] * moments[0];
                auSum += a[k] * moments[1];
                avSum += a[k] * moments[2];
            }

            double tm = auSum / acSum;
            double tv = avSum / acSum - tm * tm;
            if (tv > 0 && !double.IsInfinity(tv) && !double.IsNaN(tm))
            {
                mean = tm;
                variance = tv;
            }
        }

        public double getProbability()
        {
            return GaussianIntegral.NormCdf(mean / Math.Sqrt(1 + variance));
        }
    }

    public static class GaussianIntegral
    {
        #region Static Methods
        /// <summary>
        /// Computes the complementary error function. This function is defined 
        /// by 2/sqrt(pi) * integral from x to infinity of exp (-t^2) dt.
        /// </summary>
        /// <param name="x">The upper integration boundary.</param>
        /// <returns>2/sqrt(pi) * integral from x to infinity of exp (-t^2) dt</returns>
        private static double ErfC(double x)
        {
            if (Double.IsNegativeInfinity(x))
                return (2.0);
            else if (Double.IsPositiveInfinity(x))
                return (0.0);
            else
            {
                double z = Math.Abs(x);
                double t = 1.0 / (1.0 + 0.5 * z);
                double res = t * Math.Exp(-z * z - 1.26551223 +
                    t * (1.00002368 +
                    t * (0.37409196 +
                    t * (0.09678418 +
                    t * (-0.18628806 +
                    t * (0.27886807 +
                    t * (-1.13520398 +
                    t * (1.48851587 +
                    t * (-0.82215223 + t * 0.17087277)))))))));
                if (x >= 0.0) return (res); else return (2.0 - res);
            }
        }

        /// <summary>
        /// Computes the inverse of the complementary error function.
        /// </summary>
        /// <param name="y">A number between 0 and 2.</param>
        /// <returns>The inverse of the complementary error function</returns>
        private static double ErfCInv(double y)
        {
            if (y < 0.0 || y > 2.0)
                throw (new ArgumentOutOfRangeException("y", "Inverse complementary function not defined outside [0,2]."));
            else if (y == 0.0)
                return (Double.PositiveInfinity);
            else if (y == 2.0)
                return (Double.NegativeInfinity);
            else
            {
                double x = 0.0;

                if (y >= 0.0485 && y <= 1.9515)
                {
                    double q = y - 1.0;
                    double r = q * q;
                    x =
                        (((((0.01370600482778535 * r - 0.3051415712357203) * r + 1.524304069216834) * r - 3.057303267970988) * r + 2.710410832036097) * r - 0.8862269264526915) * q /
                    (((((-0.05319931523264068 * r + 0.6311946752267222) * r - 2.432796560310728) * r + 4.175081992982483) * r - 3.320170388221430) * r + 1.0);
                }
                else if (y < 0.0485)
                {
                    double q = Math.Sqrt(-2.0 * Math.Log(y / 2.0));
                    x =
                        (((((0.005504751339936943 * q + 0.2279687217114118) * q + 1.697592457770869) * q + 1.802933168781950) * q + -3.093354679843504) * q - 2.077595676404383) /
                    ((((0.007784695709041462 * q + 0.3224671290700398) * q + 2.445134137142996) * q + 3.754408661907416) * q + 1.0);
                }
                else if (y > 1.9515)
                {
                    double q = Math.Sqrt(-2.0 * Math.Log(1.0 - y / 2.0));
                    x =
                        (-(((((0.005504751339936943 * q + 0.2279687217114118) * q + 1.697592457770869) * q + 1.802933168781950) * q + -3.093354679843504) * q - 2.077595676404383) /
                     ((((0.007784695709041462 * q + 0.3224671290700398) * q + 2.445134137142996) * q + 3.754408661907416) * q + 1.0));
                }

                double u = (ErfC(x) - y) / (-2.0 / Math.Sqrt(Math.PI) * Math.Exp(-x * x));
                return (x - u / (1.0 + x * u));
            }
        }

        /// <summary>
        /// Computes the cummulative Gaussian distribution.
        /// </summary>
        /// <param name="t">Point of interest.</param>
        /// <returns>The cummulative Gaussian distribution at point of interest.</returns>
        public static double NormCdf(double t)
        {
            return ((ErfC(-t / 1.4142135623730951)) / 2.0);
        }

        /// <summary>
        /// Computes the Gaussian density.
        /// </summary>
        /// <param name="t">Point of interest.</param>
        /// <returns>Gaussian density at point of interest.</returns>
        public static double NormPdf(double t)
        {
            return (0.398942280401433 * Math.Exp(-(t * t / 2.0)));
        }

        /// <summary>
        /// Computes the inverse of the cummulative Gaussian distribution (qunatile function).
        /// </summary>
        /// <param name="p">Point of interest.</param>
        /// <returns>Quantile function at point of interest.</returns>
        public static double NormInv(double p)
        {
            return (-1.4142135623730951 * ErfCInv(2.0 * p));
        }

        /// <summary>
        /// Computes the additive correction of a single-sided truncated Gaussian 
        /// with unit variance.
        /// </summary>
        /// <param name="t">The point of interest.</param>
        /// <returns>Additive correction of a single-sided truncated Gaussian at point of interest.</returns>
        public static double V(double t)
        {
            double denom = NormCdf(t);

            if (denom < 2.222758749e-162)
                return (-t);
            else
                return (NormPdf(t) / denom);
        }

        /// <summary>
        /// Computes the multiplicative correction of a single-sided truncated Gaussian 
        /// with unit variance.
        /// </summary>
        /// <param name="t">The point of interest.</param>
        /// <returns>Multiplicative correction of a single-sided truncated Gaussian at point of interest.</returns>
        public static double W(double t)
        {
            double denom = NormCdf(t);
            if (denom < 2.222758749e-162)
                return ((t < 0.0) ? 1.0 : 0.0);
            else
            {
                double vt = V(t);
                return (vt * (vt + t));
            }
        }
        #endregion
    }

    public class Arguments
    {
        // Variables
        private StringDictionary Parameters;

        // Constructor
        public Arguments(string[] Args)
        {
            Parameters = new StringDictionary();
            Regex Spliter = new Regex(@"^-{1,2}|^/|=|:", RegexOptions.IgnoreCase | RegexOptions.Compiled);
            Regex Remover = new Regex(@"^['""]?(.*?)['""]?$", RegexOptions.IgnoreCase | RegexOptions.Compiled);
            string Parameter = null;
            string[] Parts;

            // Valid parameters forms:
            // {-,/,--}param{ ,=,:}((",')value(",'))
            // Examples: -param1 value1 --param2 /param3:"Test-:-work" /param4=happy -param5 '--=nice=--'
            foreach (string Txt in Args)
            {
                // Look for new parameters (-,/ or --) and a possible enclosed value (=,:)
                Parts = Spliter.Split(Txt, 3);
                switch (Parts.Length)
                {
                    // Found a value (for the last parameter found (space separator))
                    case 1:
                        if (Parameter != null)
                        {
                            if (!Parameters.ContainsKey(Parameter))
                            {
                                Parts[0] = Remover.Replace(Parts[0], "$1");
                                Parameters.Add(Parameter, Parts[0]);
                            }
                            Parameter = null;
                        }
                        // else Error: no parameter waiting for a value (skipped)
                        break;
                    // Found just a parameter
                    case 2:
                        // The last parameter is still waiting. With no value, set it to true.
                        if (Parameter != null)
                        {
                            if (!Parameters.ContainsKey(Parameter)) Parameters.Add(Parameter, "true");
                        }
                        Parameter = Parts[1];
                        break;
                    // Parameter with enclosed value
                    case 3:
                        // The last parameter is still waiting. With no value, set it to true.
                        if (Parameter != null)
                        {
                            if (!Parameters.ContainsKey(Parameter)) Parameters.Add(Parameter, "true");
                        }
                        Parameter = Parts[1];
                        // Remove possible enclosing characters (",')
                        if (!Parameters.ContainsKey(Parameter))
                        {
                            Parts[2] = Remover.Replace(Parts[2], "$1");
                            Parameters.Add(Parameter, Parts[2]);
                        }
                        Parameter = null;
                        break;
                }
            }
            // In case a parameter is still waiting
            if (Parameter != null)
            {
                if (!Parameters.ContainsKey(Parameter)) Parameters.Add(Parameter, "true");
            }
        }

        // Retrieve a parameter value if it exists
        public string this[string Param]
        {
            get
            {
                return (Parameters[Param]);
            }
        }
    }
}
